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Abstract 

The steady state properties of the mean density pop- 
ulation of infected cells in a viral spread is simulated by 
a general forest fire like cellular automaton model with 
two distinct populations of cells ( permissive and resis- 
tant ones) and studied in the framework of the mean 
field approximation. Stochastic dynamical ingredients 
are introduced in this model to mimic cells regenera- 
tion (with probability p) and to consider infection pro- 
cesses by other means than contiguity (with probability 
f). Simulations are carried on a L x L square lattice 
considering the eigth first neighbors. The mean den- 
sity population of infected cells {Di) is measured as 
function of the regeneration probability p, and anal- 
ized for small values of the ratio f/p and for distinct 
degrees of the cell resistance. The results obtained by 
a mean field like approach recovers the simulations re- 
sults. The role of the resistant parameter R {R > 2) 
on the steady state properties is investigated and dis- 
cussed in comparision with the R — 1 monocell case 
which corresponds to the self organized critical forest 
fire model. The fractal dimension of the dead cells ul- 
cers contours were also estimated and analised as func- 
tion of the model parameters.]^ 

I. INTRODUCTION 

In this work we study the steady state properties 
of a simple discrete dynamical stochastic cellular 
automaton model describing the virus spread. An 
earlier deterministic version of the present model 
was proposed by Landini et al ^ in attempt to 
study the morphology of Herpes Simplex Virus 
(HSV) corneal ulcers. In this latter model two 
distinct populations of epithelial cells were consid- 
ered accordingly with the degree of permissivity 
to infection by HSV. The susceptibility to viral 
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infection by contiguity were distinguished by as- 
suming that the permissive cells become infected 
by the presence of one or more infected neighbors 
while for the resistant ones at least R {R > 2) 
infected neighbors cells are required to spread de 
virus. Landini et al 0] considering a square lattice 
monolayer cellular automaton model (taking into 
account the first eight nearest neighbors) showed 
that the HSV ulcers can be quantitatively well 
characterized by the fractal dimension of their con- 
tour. Furthermore, they shed some light to un- 
derstand that the complex shape and evolution of 
the HSV ulcers may depend on intrinsic character- 
istics (susceptibility to infection) and topological 
distribution of the ephithelial cells. For instance, 
when the resistance parameter R = 5, Landini's 
model exhibit a dramatic changes from dendritic 
to amoeboid morphology when the concentration 
of permissive cells reaches the percolation thresh- 
old. 

The aim of this work is to investigate the dy- 
namical features of the HSV ulcers evolution, that 
were well described by the static model proposed 
by Landini a/ We are mainly concern with 
the overall distribution of ulcers in the steady 
state that should occurs after a primary infection. 
Nowadays it is well known that many viruses have 
evolved mechanisms to avoid the immune system, 
like e.g. adenovirus, murine and HSV. For instance 
Hill et al ^ describe a new mechanism by which 
HSV may evade the immune control. For those 
viruses after a primary infection the virus enters 
a latent phase and can reactivate later on at any 
time due to other factors. Therefore we introduce 
stochastic dynamical ingredients to the model pro- 
posed by Landini et al 0] to mimic regeneration of 
the ephithelial tissue and to consider re-infection 
processes due to latent virus phase. With these in- 
gredients our model can also describe a forest fire 
phenomena with two distinct species of trees with 
different degrees of resistance to burning. Actu- 
ally, the particular case with R — 1 (one specie 
of tree) recovers the Self-Organized Critical forest 
fire model with lightning probability introduced by 
Drossel and Schwabl ||]. 

In the present work, two methods are employed 
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to study the role of the resistance parameter R on 
the steady state properties of the dead ceU popula- 
tion (destroyed trees or empty sites in the counter- 
part forest fire model): (a) an analytical mean field 
approach and (b) a numerical simulation. In sec- 
tion II we describe our cellular automaton model 
and develops the mean field approach. Section III 
is devoted to discuss the numeric simulation pro- 
cedure. The results and discussions are given on 
section IV. 



II. THE MODEL AND THE MEAN FIELD 
APPROACH 

Our stochastic discrete cellular automaton is de- 
fined on a square lattice with sites, the sites be- 
ing occupied by a permissive cell, a resistant cell, 
an infected cell or a dead cell. The system is par- 
allel updated and governed by the following rules 
during one time step: 

(a) a permissive cell becomes infected by conti- 

guity if exist one or more infected neighbor 
cells on it environment. 

(b) a resistant cell becomes infected by contiguity 

if exist R or more infected neighbors cells on 
it environment. 

(c) an infected cell becomes a dead cell. 

(d) a living cell (permissive or resistant) may be- 

comes infected with probability / if there are 
no enough infected neighbors on it environ- 
ment. 

(e) a dead cell may regenerates with probability p, 

a fraction q being permissive cells and (1-?) 
being resistant cells. 

(f) the cell environment is defined by its eight 

nearest neighbors. Therefore R ranges from 
2 to 8. 

To study the mean field behavior of the sys- 
tem we will assume that the system size L is large 
enough to prevent finite-size effects, and that the 
system reaches a steady state after a transient pe- 
riod from arbitrary initial conditions. We also as- 
sume that this steady state depends only on the 
model parameters and is characterized by mean 
values of the densities of the permissive cells Dp, 
of the resistant cells Dr, of the infected cells Di 
and that of the dead cells Dd- These densities are 
constrained by normalization condition: 



Dp + Dr + D, + Dd = l 



(1) 



Now let's consider the variation of these densi- 
ties within one time step. The rate equation for 
dead cells is clearly given by 

/^Dd - A - pDd, (2) 

while for the permissive cells we have 

Ai?p = qpDd - [6 + /(I - Q)]Dp (3) 

In Eq. ^, the first term gives the mean number 
of new regenerate permissive cells and the second 
reads for the permissive cells becoming infected on 
the next time step (rules (a) and (d)), 8 = 1 — (1 — 
Di)^ being the probability for a given permissive 
cell to have one or more infected neighbors. By 
analogy for the rate of resistant cells we have 

^Dr = (1 - q)pDd - [^R + /(I - $fl)]i?r (4) 

where $ij is the probability that a given resistant 
cell is surrounded by R or more infected cells, that 
is 



n=R. 



(5) 



Finally the rate of the infected cells is given by, 
AA - (OA + A) + 



[(i-e)A + (i-*i^)-Dr]/- A 



(6) 



In the above equation the first term describes the 
number of cells infected by contiguity in one time 
step, the second gives the ones infected by other 
means and the third clearly stands for the number 
of infected cells that will become dead in the next 
time step. 

The steady state is characterized by constant 
mean densities of all populations of cells, that is, 

A A = A A = AL»j = ADrf = (7) 

In that conditions we obtain from Eqs. (0-0) that, 



Dr = 



Dd = -A 
P 



qD^ 



e + /(i-e) 



(i-g) A 



(8) 



(9) 



(10) 
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Now substituting these expressions onto the nor- 
mahzation conditions (Eq.(^) we finally obtain an 
equation for the infected cells populations as func- 
tion of the model parameters (provide f,p ^ 0) 
given by, 



1 

V 



(l-/)e + / (1-/)$^ + / 



= 1 



(11) 



This equation, which express the mean field be- 
havior of the system, can be solved numerically, 
and the figures for Di can be used to obtain the 
others steady state populations given by Eqs. (||- 
0). 



III. NUMERICAL SIMULATION 
PROCEDURE. 

We carry on numerical simulations of the present 
model by considering an L x L square lattice with 
periodic boundary conditions. The initial state is 
prepared by randomly distributing permissive cells 
with probability 0.4 (arbitrary) and resistant cells 
with probability 0.6. An infection cell is triggered 
on central site of the lattice. After a transient time 
interval (time steps) of the order of L we mea- 
sured all kind of cells densities and then average 
these values over the next L consecutive configu- 
rations. During all processes infected cells were 
triggered randomly at 1// time steps and regener- 
ated cells were randomly allowed at 1 j-p time steps. 
As far as we are interested to investigate the possi- 
bility of occurrence of self-organized critical steady 
states we fixed the model parameters such that 
1/ f ^ 1/p L . This is a necessary condition 
to allow that all living cells of an infected cluster 
become infected before new regenerated cells ap- 
pears on the clusters ends. Note that L/2 is the 
maximum time interval for infection to spread in 
a large living cells cluster. The fractal dimension 
of the contours of dead cells clusters were also es- 
timated by the box counting method and averaged 
over all simulations. 



IV. RESULTS AND DISCUSSION 

Since we are mainly interested in the properties 
of the dynamic steady state of the virus spread 
we direct our attention to the mean density of the 
population of infected cells Di or mean density of 
fire in the counterpart forest fire model. In figure 



1, we show the dependence of Di as function of 
the regeneration probability p obtained by simula- 
tion of the model on a L = 100 square lattice with 
the degree of resistance varying from R = 1 to 8, 
for the ratio f/p = 0.1 and q = 0.6. We notice 
that the curve for i? = 2 looks qualitatively simi- 
lar to the one obtained for R = 1 {SOC forest fire 
model) and is very distinct to the ones obtained 
for i? > 4. Those latter curves show an universal 
like dependence for lower values of p. Furthermore, 
the plot for i? = 3 exhibit a crossover behavior be- 
tween the R < 2 and R > A curves as p varies 
from lower to upper values. Now we make com- 
parison with the mean field like results obtained 
from Eq. 11. In figure 2 we display the results ob- 



tained by the mean field procedure recovering the 
ones obtained by simulations. For i? = 1 and 2 the 
mean field curves are in qualitative agreement of 
the one obtained by simulations, while for R = 3 
one also observes the cross-over between R < 2 
and i? > 4 behaviors with a very sharp changes 
at small but definite value of p. For i? > 4 the 
agreement between both approaches is quantita- 
tively established. 

In figures 3 and 4 we present the behavior of 
the density of dead cells (ulcers) as function of the 
regeneration probability as obtained by simulation 
and mean field approaches, respectively. The same 
qualitative and quantitative behavior observed for 
the density of infected cells (figures 1 and 2) occurs 
for the dead cell density. To explore the crossover 
behavior observed for the R — 3 case, we plotted in 
figure 5 the density of infected cells Di against the 
regeneration probability p for several values of the 
parameter if Ip) obtained by the mean field ap- 
proach (Eq. ^I]). We observe that the mean field 
solution (that is independent of finite size effects) 
changes dramatically as p approaches from a cer- 
tain value which is f/p dependent. As p decreases 
and approaches to zero the condition f jp ^ 1 is no 
longer fulfilled and the SOC state disappears 
This means that clusters of resistant cells remains 
alive for long time intervals stopping the virus 
spread. In figure 6, we compare R — Z case ob- 
tained by both methods by showing Di x p plot 
for small values of p. This plot indicates that for 
i? = 3 and for very small values of p and for a fixed 
ratio f/p the mean field predicted behavior is con- 
firmed by simulations, even for finite small lattices. 
Therefore the expected SOC state for the dynam- 
ical behavior of the clusters of dead cells (ulcers), 
similar to the one observed in the forest fire model 
(i? = 1) 1^ for small values of p should occurs only 
for R = 2 case. 
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We have also investigated the dependence of the 
mean density of infected cells Di against the varia- 
tion of the parameter q (the fraction of permissive 
regenerate cells) for a fixed value of p. This is 
shown in figure 7 from a simulation on a L = 200 
square lattice and for f/p = 0.1 and p — 0.6 . One 
observe the crossover behavior for the i? = 3 case 
as q increases from lower to upper values. 

Finally we have investigated the dependence of 
the fractal dimension Dp of the countour of the 
dead cells clusters as a function of the model pa- 
rameters. The clusters were generated by simu- 
lations and after the transient time interval ( of 
order of L ) the fractal dimension of ulcers con- 
tours were calculated by using the box counting 
method and averaged over the next L equilibrium 
configuration. 

In figure 8 we show Dp as function of small val- 
ues of (f/p) obtained for p = 0.1 and for R = 2,3 
and 4. The figures for i? > 4 as quite similar to 
the ones for i? = 4. We notice that for the models 
with resistant cells with R > 3 the fractal dimen- 
sion of the ulcers are very sensitive to small values 
of the {f/p) parameter. Under these conditions 
the densities of dead cells is of order of 0.2 (see 
figure 3) indicating that a steady state regime is 
achieved where small ulcers of very rough contour 
remain dynamically stable. For {f/p) close to 0.01 
the ulcers frontiers are less rough for i? > 2 with 
Dp ~ 1.4. This figure is close to the one estimated 
by box counting digitalized ulcers images of den- 
tritic small real ulcers (Feret's diameter ranging 
from 1.6 to 3.2 mm) 0| . On the other hand if the 
cell resistance is low {R < 2) this latter picture 
does not hold and the ulcers is allowed to spread 
over the whole tissue with high density and with 
Dp - 2. 

In conclusion, we show that our present dynam- 
ical model indicates that the degree of resistance 
of the cell to infection by contiguity should plays 
an important role on the ulcer evolution. 
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FIG. 1. Mean density of infected cells Di obtained 
by simulation for a L = 100 square lattice as function 
of the regeneration probability p. R=l — 8,f/p = 0.1, 
q — 0.6. Dashed lines are guide for the eyes. 
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FIG. 2. Mean density of infected cells Di obtained 
by mean field approximation, as function of the regen- 
eration probability p. Same parameters of figure 1. 
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FIG. 3. Mean density of dead eells Dd obtained by 
simulation on a L = 100 square lattice as function of 
the regeneration probability p. Same parameters of 
figure 1. 
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FIG. 5. Mean density of infected cells Di obtained 
by the mean field approach as function of the regener- 
ation probability p for = 3 for several values from 
the ratio f/p = 0.0001 to 1.0, with q = 0.6. 
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FIG. 4. Mean density of dead cells obtained by mean 

field approximation, as function of the regeneration 
probability. Same parameters of figure 1. 
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FIG. 6. Mean density of infected cells Di as func- 
tion of the regeneration probability p for resistance 
parameter R = 3. Full line is obtained from mean 
field approach and open circles indicates figures from 
simulations. Same parameters of figure 1. 



5 



Simulation 



References 



.5 0.2 

S 



0- 


-GR = 2 




-□R = 3 


0- 


-OR = 4 



^0 

0'' 



. -a' 



[1] G. Landini, G.P. Mission and P.I. Murray, Fractals 
in the Na.tural and Applied Sciences, M.M. Novak 
Ed., Elsevier Science B.V., North Holland (1994); 

[2] A. Hill, P. Jugovic, I. York, G. Russ, J. Bennink, J. 
Yewdell, H. Ploegh and D. Johnson, Nature, 375, 
411, (1995); 

[3] B. Drossol and F. Schwabl, Phys. Rev. Lett. 69, 
1629, (1992); 

[4] B. Drossel, S. Clar and F. Schwabl, Phys. Rev. E 
50, 1009, (1994). 



0.0 0.2 0.4 0.6 0.8 1.0 
Fraction of permissive regenerate cells 

FIG. 7. Mean density of infected cells Di obtained 
by simulation on a L = 200 square lattice as function 
of the fraction of regeneration of permissive cells q for 
f/p = 0.1 and p = 0.6. 
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FIG. 8. Fractal dimension of ulcers contours Dp ob- 
tained by simulation as function of the f/p ratio for 
p = 0.01 and q = 0.6 on a L = 100 square lattice, for 
i? = 2, 3 and 4. 
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